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ABSTRACT 

We present a series of high-resolution (20-2000 Mq, 0.1-4 pc) cosmological zoom-in simula¬ 
tions at z > 6 from the Feedback In Realistic Environment (FIRE) project. These simulations 
cover halo masses 10^-10^^ Mq and rest-frame ultraviolet magnitude Muv = —9 to —19. 
These simulations include explicit models of the multi-phase ISM, star formation, and stellar 
feedback, which produce reasonable galaxy properties at z = 0-6. We post-process the snap¬ 
shots with a radiative transfer code to evaluate the escape fraction (/esc) of hydrogen ionizing 
photons. We find that the instantaneous /esc has large time variability (0.01%-20%), while the 
time-averaged /esc over long time-scales generally remains ^5%, considerably lower than the 
estimate in many reionization models. We find no strong dependence of /esc on galaxy mass 
or redshift. In our simulations, the intrinsic ionizing photon budgets are dominated by stellar 
populations younger than 3 Myr, which tend to be buried in dense birth clouds. The escap¬ 
ing photons mostly come from populations between 3-10 Myr, whose birth clouds have been 
largely cleared by stellar feedback. However, these populations only contribute a small frac¬ 
tion of intrinsic ionizing photon budgets according to standard stellar population models. We 
show that /esc can be boosted to high values, if stellar populations older than 3 Myr produce 
more ionizing photons than standard stellar population models (as motivated by, e.g., models 
including binaries). By contrast, runaway stars with velocities suggested by observations can 
enhance /esc by only a small fraction. We show that “sub-grid” star formation models, which 
do not explicitly resolve star formation in dense clouds with \ cm“^, will dramatically 
over-predict /esc- 

Key words: galaxies: formation - galaxies: evolution - galaxies: high-redshift - cosmology: 
theory 


1 INTRODUCTION 

Star-forming galaxies at high redshifts are thought to be the dom¬ 
inant source of hydrogen reionization (e.g. Madau et al. 1999; 
Faucher-Giguere et al. 2008; Haardt & Madau 2012). Therefore, 
the escape fraction of hydrogen ionizing photons (/esc) from these 
galaxies is an important, yet poorly constrained, parameter in un¬ 
derstanding the reionization history. 

Models of cosmic reionization are usually derived from the 
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galaxy ultraviolet luminosity function (UVLF; e.g. Bouwens et al. 
2011; McLure et al. 2013), Thomson scattering optical depths in¬ 
ferred from Cosmic Microwave Background (CMB) measurements 
(Hinshaw et al. 2013; Planck Collaboration et al. 2013), Lya for¬ 
est transmission (e.g. Fan et al. 2006). They often require high /esc 
in order to match the ionization state of the intergalactic medium 
(IGM) by z = 6 (e.g. Ouchi et al. 2009; Kuhlen & Faucher-Giguere 
2012; Finkelstein et al. 2012; Robertson et al. 2013). For exam¬ 
ple, Finkelstein et al. (2012) and Robertson et al. (2013) suggested 
/esc >13% and /esc > 20%, respectively, assuming all the ioniza¬ 
tion photons are contributed by galaxies brighter than Muv = —13. 
However, such constraints on /esc are always entangled with the 
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uncertainties at the faint end of UVLF, since low-mass galaxies 
can play a dominant role in providing ionizing photons due to their 
dramatically increasing numbers. For example, Finkelstein et al. 
(2012) derived that reionization requires a much higher escape frac¬ 
tion /esc > 34% if one only accounts for the contribution of galax¬ 
ies brighter than Muv = —18. Also, Kuhlen & Faucher-Giguere 
(2012) showed that even applying a cut off on UV magnitude at 
Muv = —13, the required escape fraction at z = 6 varies from 6%- 
30% when changing the faint-end slope of UVLF within observa¬ 
tional uneertainties. Furthermore, it is also not clear how /esc de¬ 
pends on galaxy mass and evolves with redshift, which makes the 
problem more complicated. 

Therefore, independent constraints on /esc are necessary to 
disentangle these degeneraeies. Star-forming galaxies at lower red- 
shifts should provide important insights into their high-redshift 
counterparts. In the literature, high escape fractions from 10% up to 
unity have been reported in various samples of Lyman break galax¬ 
ies (LBGs) and Lya emitters (LAEs) around z ~ 3 (e.g. Steidel et 
al. 2001; Shapley et al. 2006; Vanzella et al. 2012; Nestor et al. 
2013). These measurements are based on detection of rest-frame 
Lyman continuum (LyC) emission from either individual galaxies 
or stacked samples, so the exact value of /esc depends on uncertain 
dust and IGM attenuation correction. Similar observations at lower 
redshifts always show surprisingly low escape fractions. In the lo¬ 
cal universe, the only two galaxies which have confirmed LyC de¬ 
tection suggest /esc to be only ~ 2%-3% (Leitet et al. 2011, 2013). 
At z ~ 1, stacked samples have been used to derive upper limits 
as low as /esc < l%-2% (e.g. Cowie et al. 2009; Siana et al. 2010; 
Bridge et al. 2010). Even at z ~ 3, low escape fractions (< 5%) 
have also been reported in some galaxy samples (e.g. Iwata et al. 
2009; Boutsia et al. 2011). Recent careful studies have revealed that 
a considerable fraction of specious LyC detection at z ~ 3 is due to 
contamination from foreground sources (Vanzella et al. 2010; for a 
very recent study, see Siana et al. 2015), which could at least partly 
account for the apparent contradiction between these observations. 
Nevertheless, given the large uneertainty in these studies, no con¬ 
vincing conclusion can be reached so far from current observations. 

Previous numerical simulations of galaxy formation also pre¬ 
dict a broad range of /esc, and even contradictory trends of the de¬ 
pendence of /esc on halo mass and redshift. Eor example, Razoumov 
& Sommer-Larsen (2010) found /esc decreases from unity to a few 
precent with increasing halo mass from 10^ ^-10^^ ^ M©. Similarly, 
Yajima et al. (2011) also found their /esc decreases from 40% at 
halo mass 10^ M© to 7% at halo mass 10^^ M©. On the other 
hand, Gnedin et al. (2008) found increasing /esc with halo mass 
in lO^^-lO^^ M©. They also reported significantly lower escape 
fraction of 1% — 3% for the most massive galaxies in their sim¬ 
ulations and <0.1% for the smaller ones. Razoumov & Sommer- 
Larsen (2010) also found /esc decreases from z = 4-10 at fixed halo 
mass, while Yajima et al. (2011) found no dependence of /esc on 
redshift. At lower masses. Wise & Cen (2009) found /esc ~ 5%- 
40% and /esc ~ 25%-80% by invoking a normal initial mass func¬ 
tion (IME) and a top-heavy IME, respectively, for galaxies of halo 
mass in lO^'^-lO^'^ M©; whereas Paardekooper et al. (2011) re¬ 
ported lower escape fraction of 10“^-0.1 in idealized simulations 
of galaxy masses 10^-10^ M©. 

Most of the intrinsic ionizing photons are produced by mas¬ 
sive stars of masses in 10-100 M©, which are originally born in 
giant molecular clouds (GMCs). The majority of the ionizing pho¬ 
tons are instantaneously absorbed by the dense gas in the GMCs 
and generate H II regions. These “birth clouds” must be disrupted 
and dispersed by radiation pressure, photoionization, H II thermal 


pressure, and supernovae before a eonsiderable fraetion of ionizing 
photons are able to escape (e.g. Murray et al. 2010; Kim et al. 2013; 
Paardekooper et al. 2015). Therefore, to study the escape fraction 
of ionizing photons using simulations, one must resolve the multi¬ 
phase structure of the interstellar medium (ISM) and “correctly” 
describe star formation and stellar feedback. Many previous sim¬ 
ulations adopt very approximate or “sub-grid” ISM and feedback 
model, which can lead to many differenees between those stud¬ 
ies. Recent studies have noted the importance of resolving the ISM 
structure around the stars and started to adopt more detailed treat¬ 
ments of the ISM and stellar feedback physics (Kim et al. 2013; 
Kimm & Cen 2014; Wise et al. 2014; Paardekooper et al. 2013, 
2015). Eor example. Wise et al. (2014) performed radiative hydro- 
dynamical simulations with state-of-art ISM physics and chemistry, 
star formation, and stellar feedback models and found /esc drops 
from 50% to 5% with increasing halo mass in 10^-10^ ^ M© at 
z > 7. They conclude that more massive galaxies are not likely to 
have high eseape fractions, but are unable to simulate more massive 
systems. Kimm & Cen (2014) explored more physieally motivated 
models of supernovae (SN) feedback and found average escape 
fraction of ~ 11% for galaxies in 10^-10^^'^ M©. Paardekooper et 
al. (2015) argued that the dense gas within 10 pc from young stars 
provides the main constraint on the escape fraction. They found 
in their simulation that about 70% of the galaxies of halo mass 
above 10^ M© have escape fraction below 1%. But none of these 
simulations has been run to z = 0 to confirm that the models for 
star formation, feedback, and the ISM produee reasonable results 
in eomparison to observations. 

The Eeedback in Realistic Environment (EIRE) project^ (Hop¬ 
kins et al. 2014) is a series of cosmological zoom-in simulations 
that are able to follow galaxy merger histories, interaction of galax¬ 
ies with IGM, and many other processes. The simulations include 
a full set of realistic models of the multi-phase ISM, star forma¬ 
tion, and stellar feedbaek. The first series of EIRE simulations 
run down to z = 0 reproduce reasonable star formation histories, 
the stellar mass-halo mass relation, the Kennicutt-Schmidt law, 
and the star-forming main sequence, for a broad range of galaxy 
masses (M* = lO'^-lO^^ M©) from z = 0-6 (Hopkins et al. 2014). 
Cosmological simulations with the EIRE stellar feedback physics 
self-consistently generate galaetic winds with velocities and mass 
loading factors broadly consistent with observational requirements 
(Muratov et al. 2015) and are in good agreement with the observed 
covering fractions of neutral hydrogen in the halos of z = 2-3 
LBGs (Eaucher-Giguere et al. 2014). In previous studies of iso¬ 
lated galaxy simulations, these models have also been shown to 
reproduce many small scale observations, including the observed 
multi-phase ISM structure, density distribution of GMCs, GMC 
lifetimes and star formation efficiencies, and the observed Larson’s 
law sealings between cloud sizes and struetural properties, from 
scales < 1 pc to > kpc (e.g. Hopkins et al. 2011, 2012). A realistic 
model with these properties is necessary to study the production 
and propagation of ionizing photons inside a galaxy. 

In this work, we present a separate set of cosmological sim¬ 
ulations at z > 6, performed with the same method and models at 
extremely high resolution (particle masses 20-2000 M©, smooth¬ 
ing lengths 0.1-4 pc). These simulations cover galaxy halo masses 
10^-10^^ M© and rest-frame ultraviolet magnitudes Muv = —9 to 
—19 at z = 6. We then evaluate the escape fraction of ionizing pho¬ 
tons with Monte Carlo radiative transfer calculations. We describe 


^ FIRE project website: http://fire.northwestern.edu 
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the simulations and present the properties of our galaxies in Sec¬ 
tion 2 and 3. In Section 4, we describe the Monte Carlo radiative 
transfer code and compile our main results on the escape fractions 
and there dependence on galaxy mass and cosmic time. In Section 
5, we show how the UV background and star formation prescrip¬ 
tions affect our results. We also discuss the effects of runaway stars 
and extra ionizing photon budgets contributed by intermediate-age 
stellar populations, as motivated by recent observations and stellar 
models. We summarize and conclude in Section 6 . 


2 THE SIMULATIONS 

This work is part of the FIRE project (Hopkins et al. 2014). All the 
simulations use the newly developed GIZMO code (Hopkins 2014) 
in “P-SPH” mode. P-SPH adopts a Lagrangian pressure-entropy 
formulation of the smoothed particle hydrodynamics (SPH) equa¬ 
tions (Hopkins 2013), which eliminates the major differences be¬ 
tween SPH, moving-mesh, and grid codes, and resolves many well- 
known issues in traditional density-based SPH formulations. The 
gravity solver is a heavily modified version of the GADGET-3 
code (Springel 2005); and P-SPH also includes substantial im¬ 
provements in the artificial viscosity, entropy diffusion, adaptive 
time-stepping, smoothing kernel, and gravitational softening algo¬ 
rithm. We refer to Hopkins (2013, 2014) for more details on the 
numerical recipes and extensive test problems. A list of the simu¬ 
lations in this work is presented in Table 1, while the parameters 
there will be introduced in the rest of this section. 

The simulations in this work are of a separate series from other 
EIRE simulations. A large cosmological box was first simulated at 
low resolution down to z = 5, and then halos of masses in 10^- 
lO^^M© at that time were picked and re-simulated in a smaller box 
at much higher resolution with the multi-scale “zoom-in” initial 
conditions generated with the MUSIC code (Hahn & Abel 2011), 
using second-order Lagrangian perturbation theory. The resolution 
of the simulations in this work can be roughly divided into two 
categories, which we refer as “high resolution” (HR) and “medium 
resolution” (MR), although the specific initial particle mass may 
vary according to the size of the system. Some initial conditions 
we adopt in the simulations and general properties of the galaxies 
at z = 6 are listed in Table 1. We will show in Section 3 that they 
are typical in most of their properties and thus can be considered as 
“representative” in this mass range. 

In our simulations, gas follows an ionized-i-atomic-i-molecular 
cooling curve from 10-10^^ K, including metallicity-dependent 
fine-structure and molecular cooling at low temperatures and high- 
temperature metal-line cooling followed species-by-species for 11 
separately tracked species (Wiersma et al. 2009a). We do not in¬ 
clude a primordial chemistry network nor try to model the forma¬ 
tion of Pop III stars, but apply a metallicity fioor of Z = 10 -“z© 
in the simulations. Therefore, we will focus our analysis at z < 11, 
when our galaxies are sufficiently metal-enriched. 

At each timestep, the ionization states are determined from 
the photoionization equilibrium equations described in Katz et al. 
(1996) and the cooling rates are calculated from a compilation 
of CLOUDY runs, by applying a uniform but redshift-dependent 
photo-ionizing background tabulated in Eaucher-Giguere et al. 
(2009)^, and photo-ionizing and photo-electric heating from lo- 

^ The photo-ionizing background stars to kick in at z = 10.6 and is avail¬ 
able at http://galaxies.northwestern.edu/uvb/. 


cal sources. Gas self-shielding is accounted for with a local Jeans- 
length approximation, which is consistent with the radiative trans¬ 
fer calculation in Eaucher-Giguere et al. (2010). In this work, we 
also post-process the simulations with full radiative transfer calcu¬ 
lation and re-compute the ionization states. We find consistent re¬ 
sults between the post-processing and on-the-fiy calculations (see 
Appendix A for details). 

The models of star formation (SE) and stellar feedback imple¬ 
mented in the EIRE simulations are developed and presented in a 
series of papers (Hopkins et al. 2011, 2012, 2013, 2014, and refer¬ 
ences therein). We briefiy summarize their main features here and 
refer to the references for more details and discussion. We follow 
the SE criteria developed in Hopkins et al. (2013) and allow stars to 
form only in molecular and self-gravitating gas clouds with num¬ 
ber density above some threshold Wth- We choose nth = 100 cm“^ 
as the fiducial value. It corresponds to the typical density of GMCs 
and is much larger than the mean ISM density in our simula¬ 
tions^. In zSmlOh, we set nth = 1000 cm“^ for a convergence 
test. SE occurs at 100% efficiency per free-fall time when the gas 
meets these criteria (i.e. p* = p/to). This SE prescription adap¬ 
tively selects the largest over-densities and automatically predicts 
clustered SE (Hopkins et al. 2013). It is also motivated by much 
higher-resolution, direct simulations of dense, star-forming clouds 
(Padoan & Nordlund 2011; Vazquez-Semadeni et al. 2011; Eed- 
errath et al. 2011). A star particle inherits the metallicity of each 
tracked species from its parent gas particle, and its age is deter¬ 
mined by its formation time in subsequent timesteps. 

The zSmlOe run is intentionally designed to mimic “sub-grid” 
star formation models as commonly adopted in low-resolution sim¬ 
ulations that cannot capture the star formation in dense gas clouds. 
In this run, we lower nth to 1 cm“^ and allow stars to form at 2% 
efficiency per free-fall time in all gas above 1 cm“^ but not self- 
gravitating (still 100% efficiency in self-gravitating gas). This will 
result in a wide spatial and density distribution of SE and means 
that stars do not need to form in high-density structures. 

Every star particle is treated as a single stellar population with 
known age, metallicity, and mass. Then all the quantities associated 
with feedback, including ionizing photon budgets, luminosities, 
stellar spectra, supernovae (SNe) rates, mechanical luminosities of 
stellar winds, metal yields, etc., are directly tabulated from the stel¬ 
lar population models in STARBURST99 (Leitherer et al. 1999), 
assuming a Kroupa (2002) initial mass function (IME) from 0.1- 
100 M 0 . In principle, this “IME-averaged” approximation breaks 
down in our HR simulations, where the mass of a star particle is 
only 10-100 Mq . Previous studies showed that it has little effect on 
global galaxy properties (e.g. Hopkins et al. 2014, and references 
therein). We also test and confirm that this approximation does not 
affect our results on escape fraction (see Section 4). 

We account for different mechanisms of stellar feedback, in¬ 
cluding: ( 1 ) local and long-range momentum flux from radiative 
pressure; ( 2 ) energy, momentum, mass and metal injection from 
SNe and stellar winds; and (3) photoionization and photoelectric 
heating. We apply the Type-II SNe rates from STARBURST99 and 
Type-la SNe rates following Mannucci et al. (2006), when a star 
particle is older than 3 Myr and 40 Myr, respectively. We assume 
that every SN ejecta has an initial kinetic energy of 10^^ ergs, which 
is coupled to the gas as either thermal energy or momentum, de¬ 
pending on whether the cooling radius can be resolved (see Hop- 


^ On the other hand, the threshold is much less than the highest density 
these simulations can resolve, to save computational expense. 
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Figure 1. Gas and stars in z5m09 (left column), zSmlOmr (middle column), and z5mll (right column), at z = 9.6 (upper panels) and z = 6.0 (lower panels), 
respectively. Gas images show log-weighted projected gas density. Magenta shows cold molecular/atomic gas {T < 1000 K), green shows warm ionized 
gas (10^ ^ ^ 10^ K), and red shows hot gas (T > 10^ K) (see Hopkins et al. 2014 for details). Stellar images are mock ulgjr composites. We use 

STARBURST99 to determine the SED of each star particle from its known age and metallicity, and then ray-tracing the line-of-sight flux, attenuating with a 
MW-like reddening curve with constant dust-to-metals ratio for the abundance at each point. White circles show the position and halo virial radii of each main 
galaxy (see text) identifled by the AHF code. Gas and star images of the same snapshot use the same projection and the same box size along each direction. We 
can clearly see a complicated, multi-phase ISM structure, with inflows, outflows, mergers, and star formation in dense clouds all occurring at the same time. 
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Table 1. The simulations. 


Name 

Boxsize 
(ft-‘ Mpe) 

Myir 

(Mo) 

Mb 

(Mo) 

(pc) 

^dm 

(Mo) 

Cdm 

(pc) 

nth 

(cm“^) 

Mq) 

Muv 
(AB mag) 

Resolution 

z5m09 

1 

7.6e8 

16.8 

0.14 

81.9 

5.6 

100 

3.1e5 

-10.1 

HR 

zSmlO 

4 

1.3el0 

131.6 

0.4 

655.6 

7 

100 

2.7e7 

-14.8 

HR 

zSmlOmr 

4 

1.5el0 

l.le3 

1.9 

5.2e3 

14 

100 

5.0e7 

-17.5 

MR 

zSmlOe^ 

4 

1.3el0 

l.le3 

1.9 

5.2e3 

14 

1 

2.4e7 

-16.1 

MR 

zSmlOh 

4 

1.3el0 

l.le3 

1.9 

5.2e3 

14 

1000 

6.6e7 

-16.4 

MR 

zSmll 

10 

5.6el0 

2.1e3 

4.2 

1.0e4 

14 

100 

2.0e8 

-18.5 

MR 


Initial conditions of our simulations and simulated galaxy properties at z = 6: 

(1) Name: Simulation designation. 

(2) Boxsize: Zoom-in box size of our simulations. 

(3) Mviri Halo mass of the primary galaxy at z = 6. 

(4) m}j\ Initial baryonic (gas and star) particle mass in the high-resolution region. 

(5) Minimum baryonic force softening (minimum SPH smoothing lengths are comparable or smaller). Force softening is 
adaptive (mass resolution is fixed). 

(6) rndm- Dark matter particle mass in the high-resolution regions. 

(7) Cdm- Minimum dark matter force softening (fixed in physical units at all redshifts). 

(8) «th- Density threshold of star formation (see Section 2). 

(9) M*: Stellar mass of the primary galaxy at z = 6. 

(10) Muv- Galaxy UV magnitude (absolute AB magnitude at 1500 A). 

(11) Resolution: Whether a simulation is of ultra-high resolution (HR) or of medium resolution (MR). 

Note: 

^ This simulation is intentionally designed to mimic “sub-grid” star formation models that are usually adopted in low-resolution 
cosmological simulations. We not only lower the star formation density threshold to nth = 1 cm“^, but also allow star formation 
at 2% efficiency per free-fall time if gas reaches nth but is not self-gravitating. 


kins et al. 2014; Martizzi et al. 2015, for more details). We also 
follow Wiersma et al. (2009b) and adopt Type-II SNe yields from 
Woosley & Weaver (1995) and Type-la yields from Iwamoto et al. 
(1999). We do not model SNe and metal enrichment from Pop III 
stars. 

We emphasize that the on-the-fly photoionization is treated in 
an approximate way in our simulations - we move radially out¬ 
wards from the star and ionize each nearest neutral gas particle 
until the photon budget is exhausted. This treatment allows ion¬ 
izing regions to overlap and expand, and is qualitatively reasonable 
in intense star-forming regions. However, when the gas distribu¬ 
tion is highly asymmetric around an isolated star particle, their ion¬ 
ization states might not be accurately captured in the simulations. 
Nonetheless, as we will post-process our simulations with full ra¬ 
diative transfer code to trace the propagation of ionizing photons 
and re-compute the ionization states (Section 4), this approxima¬ 
tion will have little effect on the escape fraction we evaluate. Also, 
in the region where the gas density is extremely high, photoion¬ 
ization may not be well-captured due to resolution limits. But we 
confirm that this neither has strong dynamical effect on gas struc¬ 
ture in high-density regions nor affects the escape fraction"^. 


^ In our simulations, star particles have similar mass to gas particles. Ap¬ 
plying a (Kroupa 2002) IMF and standard stellar population model, in re¬ 
gions with nth ~ 100 cm“^, the ionizing photons emitted from a young star 
particle can ionize the mass of two gas particles. However, some clouds 
reach densities > 2000 cm“^, where one needs to collect the ionizing pho¬ 
ton budgets from 10 young star particles to fully ionize a single gas particle. 
In the code, the on-the-fly estimate of HII photoionization feedback treats 
this limit stochastically (see Hopkins et al. 2011), so we might risk under¬ 
estimating the dynamical effects of photo-heating. Therefore, we run a sim¬ 
ulation where we artificially boost the ionizing photon budget by a factor of 
10, which is not physical but dramatically reduces the stochastic variations. 
We find that the typical gas density of star-forming clouds and the average 
escape fractions (computed from our post-processing radiative transfer, see 


The simulations described in Table 1 adopt a standard 
fiat ACDM cosmology with cosmological parameters Hq — 
70.2 km s“‘ Mpc"‘, S^a = 0.728, = 1 - S2a = 0.272, = 

0.0455, erg = 0.807 and n — 0.961, which are within the uncertainty 
of current observations (e.g. Hinshaw et al. 2013; Planck Collabo¬ 
ration et al. 2013). 

3 GALAXY PROPERTIES 
3.1 Halo Identification 

The galaxies in our simulations have different assembly histories at 
high redshifts. The smallest galaxy, z5m09, evolves primarily via 
accretion and passive evolution, while the more massive ones have 
undergone multiple mergers at earlier times. We use the Amiga 
Halo Finder (AHF; Gill et al. 2004; Knollmann & Knebe 2009) 
to identify halos in the simulations. The AHF code uses an adap¬ 
tive mesh refinement method. We choose the centre of a halo as 
the centre of mass of all particles in the finest refinement level 
and adopt the virial overdensity from Bryan & Norman (1998). 
In this work, we only consider the main galaxies that are well- 
resolved in the simulations. We exclude those that are contaminated 
by low-resolution paticles, not sufficiently resolved (contain less 
than 10^/10^ bound particles in MR/HR runs, or have stellar mass 
lower than 10% of the most massive galaxy in each snapshot), and 
subhalos/satellite galaxies. Some example images of gas and stars 
at different redshifts are presented in Figure 1. The white circles 
show the virial radius of each halo. As the figure shows, the more 
massive systems were assembled by merging several smaller halos 
at early time. 

Section 4) are very similar to our standard runs. Therefore, we confirm that 
the on-the-fly photoionization feedback approximation in our simulations 
does not strongly affect our results. 


© 0000 RAS, MNRAS 000, 1-17 



6 Ma et al. 


Young Middle-aged Old 



-10 12 _3 _2 -1 0 1 -6 -5 -4 -3 -2 -1 


logriH (cm lognn (cm log^n (cm 



23 2345 3 4 

logT(K) logT(K) logr(K) 


Figure 2. ISM structure in a random star neighborhood. We show density (top panels) and temperature (bottom panels) maps of a slice around a(n) young 
(~ 1 Myr, left column), middle-aged (~ 5 Myr, middle column), and old (~ 40 Myr, right column) star particle. Each box is 300 pc along each direction. The 
yellow stars represent the position of the star particle. We clearly see that young stars - which produce most of the ionizing photons - are buried in HII regions 
inside their dense birth clouds. By > 10 Myr, the clouds are totally destroyed and most sightlines to the stars have low column densities, but these stars no 
longer produce many ionizing photons. 


3.2 Multi-phase ISM Structure 

One advantage of our simulation is that we are able to explicitly re¬ 
solve a realistic multi-phase ISM, star formation, and stellar feed¬ 
back. Figure 1 shows the distribution cold, warm, and hot phase 
of gas on galactic scale. In Figure 2, we show some examples of 
ISM structure on sub-kpc scale around star particles of different 
ages from zSmlOmr. The left column is the density and tempera¬ 
ture maps around a star particle of age 1 Myr (before the first SNe 
explode at ~ 3 Myr). As expected from our star formation crite¬ 
ria, newly formed stars are embedded in their dense “birth” clouds. 
Within the central few pc around the star particle, the dense gas is 
ionized and heated by ionizing photons from the star and an H II 
region forms^. The middle column shows the ISM structure around 
an intermediate-age star particle (3-10 Myr), where there has just 
been a SN explosion (the example is 5 Myr old). The birth cloud 
has been largely dispersed and cleared by radiation pressure and SN 
feedback, opening a large covering fraction of low-density regions. 
In contrast, old star particles (right column, ~ 40 Myr) tend to be 
located in a warm, ambient medium. The ISM structures around 


^ For a typical gas density of 100 cm ^ and an ionizing photo budget 10^^-^ 
s“Mn this simulation, the Stromgren radius is around 5 pc. 


star particles of different ages are very important in understanding 
the propagation of ionizing photons. 

3.3 Galaxy Masses, Stellar Mass Assembly, and Star 
Formation History 

As has been shown in Hopkins et al. (2014), with the stellar feed¬ 
back models described here (with no tuned parameters), the FIRE 
simulations predict many observed galaxy properties from z = 0- 
6: the stellar mass-halo mass relation, the Kennicutt-Schimidt 
law, star formation histories (SFHs), and the star-forming main se¬ 
quence. The simulations in this work are of much higher resolution 
and focus on higher redshifts than those in Hopkins et al. (2014). 
We extend their analysis and present the stellar mass-halo mass 
relation at z = 6, 7, 8, and 9.6 for our simulated galaxies in Fig¬ 
ure 3. We compare our results with the simulations from Hopkins 
et al. (2014) at z = 6 and the observationally inferred relation from 
Behroozi et al. (2013) at z = 7 and 8 (note the relation in Behroozi 
et al. 2013 at z = 6 does not overlap with the halo masses presented 
here). We confirm that our simulations predict stellar masses con¬ 
sistent with observations at these redshifts. It is also reassuring that 
the stellar masses in these simulations are well converged, despite 
those from Hopkins et al. (2014) having much poorer resolution. 
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Figure 3. Galaxy stellar mass-halo mass relation at z = 6, 7, 8 , and 9.6. 
We compare the relation with the simulations from Hopkins et al. (2014) 
at z = 6 (small black dots), and the observationally inferred relation in 
Behroozi et al. (2013, z = 7.0 and z = 8.0 only, cyan lines). The black 
dotted lines represent the relation if all baryons turned into stars (i.e. 
M* = ft Mhaio)- Our simulations are broadly consistent with observations. 
These simulations are consistent with those in Hopkins et al. (2014), al¬ 
though the latter have much lower resolution. It is reassuring that the stellar 
mass is converged in these runs. 
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Figure 4. UV magnitude at 1500 A as a function of halo mass for the sim¬ 
ulated galaxies, color coded by redshift. Galaxies at z =6, 7, 8 , and 9.6 are 
shown by blue, green, red, and yellow points, respectively. The numbers 
are calculated by converting the intrinsic luminosity at 1500 A to absolute 
AB magnitude. Dotted lines show the abundance matching from Kuhlen & 
Faucher-Giguere (2012, fig 3) at z = 4 (black), 7 (magenta), and 10 (cyan). 
The simulations are qualitatively consistent with the abundance matching, 
and span the range of Muv = —9 to -19 that is believed to dominate reion¬ 
ization. 


Figure 5. Star formation rate-stellar mass relation at z = 6 , 7, 8 , and 9.6 
for the most massive galaxy in each simulation. The cyan lines illustrate the 
observed relation from a sample of galaxies at z = 6-8.7 in McLure et al. 
(2011). Our simulated galaxies agree with the observed relation where they 
connect at log(M*/M 0 ) = 8.25. 


In Figure 4, we present the relation between UV magnitude 
at rest-frame 1500 A and halo mass for our simulated galaxies at 
z = 6, 7, 8, and 9.6. To obtain the UV magnitudes, we first calculate 
the specific luminosity at 1500 A for each star particle by interpo¬ 
lating the stellar spectra tabulated from STARBURST99 as a func¬ 
tion of age and metallicity, and then convert the galaxy-integrated 
luminosity to absolute AB magnitude. In Figure 4, we also compare 
with the interpolated abundance matching from Kuhlen & Faucher- 
Gigute (2012, dotted lines). The simulations are qualitatively con¬ 
sistent with the abundance matching results, and lie within the sys¬ 
tematic observational uncertainties. Given that in this simple cal¬ 
culation, we ignore the attenuation from dust inside the galaxy and 
along the line-of-sight in the IGM, and that the abundance matching 
is very uncertain at the faint end, we do not further discuss the com¬ 
parison with these results. The simulated galaxies cover Muv = — 9 
to —19 at these redshifts, which are believed to play a dominant 
role in providing ionizing photons during reionization (e.g. Finkel- 
stein et al. 2012; Kuhlen & Faucher-Giguere 2012; Robertson et 
al. 2013). Most of these galaxies are too faint to be detectable in 
current observations; our z5mll galaxy is, however, just above the 
detection limit (Muv ~ — 18) of many deep galaxy surveys beyond 
z ~ 6. Next generation space and ground-based facilities, such as 
the James Webb Space Telescope (JWST) and the Thirty Meter Tele¬ 
scope (TMT) may push the detection limit down to Muv ~ —15.5 
(e.g. Wise et al. 2014) and many of our simulated galaxies will then 
lie above the detection limits of future deep surveys. 

Figure 5 shows the star formation rate-stellar mass relation at 
z = 6, 7, 8, and 9.6 for the most massive galaxy in each simula¬ 
tion. Our simulated galaxies agree with the observed relation from 
a sample of galaxies at z = 6-8.7 in McLure et al. (2011) where 
they connect at log(M*/Mo) = 8.25. We also present the growth 
of galaxy stellar mass and instantaneous star formation rates (SFRs) 
as a function of cosmic time for these galaxies in the top two pan¬ 
els of Figure 6 (the open symbols represent the time-averaged SFR 
on 100 Myr time-scales). All these galaxies show significant short- 
time-scale variabilities in their SFRs, associated with the dynam¬ 
ics of fountains, feedback, and individual star-forming clouds. On 
larger time-scales (e.g. 100 Myr), the fiuctuations in SFRs become 
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Table 2. Parameters used for the radiative transfer calculation. 


Name 

/min 

(pc) 

Amax 

/largest 

(pc) 

Aphoton 

Auvb 

z5m09 

25 

250 

$40 

3e7 

3e7 

z5ml0 

25 

300 

$80 

3e7 

3e7 

zSmlOmr 

50 

250 

$ 100 

3e7 

3e7 

zSmlOe 

50 

300 

$80 

3e7 

3e7 

zSmlOh 

50 

250 

$ 100 

3e7 

3e7 

z5mll 

50 

300 

$ 100 

4e7 

4e7 


( 1 ) /min' the minimum cell size. 

(2) Mnax: the maximum number of cells along each dimension. 

(3) /largest- the Cell size for the largest galaxy in the last snapshot. 

(4) A/photon- number of photon packages being transported. 


weaker and are mostly driven by mergers and global instabilities 
(see the discussion in Hopkins et al. 2014). 

It is worth noticing that our four zSmlOx simulations have 
similar global galaxy properties, despite different resolutions and 
SF prescriptions adopted in these runs. This is because the galaxy- 
averaged star formation efficiency is regulated by stellar feedback 
(1-2% per dynamical time, e.g. Hopkins et al. 2011; Ostriker & 
Shetty 2011; Agertz et al. 2013; Faucher-Giguere et al. 2013), not 
the SF criteria (Hopkins et al. 2013). However, SF criteria do affect 
the spatial and density distribution of SF. In zSmlOe, the SF density 
threshold Wth = 1 cm“^ is comparable or slightly larger than the 
mean density of the ISM, so that it can be easily reached even in the 
diffuse ISM. Also, since SF takes place at 100% efficiency per free- 
fall time once the gas becomes self-gravitating, many stars form 
just above the threshold before the gas clouds can further collapse 
to higher densities. As a consequence, stars are formed either in the 
diffuse ISM or in gas clouds of densities orders of magnitude lower 
than those in our standard runs. We emphasize that the zSmlOe run 
is not realistic but is designed to mimic star formation models as 
adopted in low-resolution simulations where the GMCs cannot be 
resolved. 


4 ESCAPE FRACTION OF IONIZING PHOTONS 

We post-process every snapshot with a three-dimensional Monte 
Carlo radiative transfer (MCRT) code to evaluate the escape frac¬ 
tion of hydrogen ionizing photons from our simulated galaxies. The 
code is derived from the MCRT code SEDONA base (Kasen et 
al. 2006) and focuses specifically on radiative transfer of hydrogen 
ionizing photons in galaxies (see Fumagalli et al. 2011, 2014). For 
each galaxy, we calculate the intrinsic ionizing photon budget for 
every star particle within Rvir to obtain the galaxy ionizing pho¬ 
ton production rate gmt- We use the Padova tracks with AGB stars 
in STARBURST99 with a metallicity Z = 0.0004 (0.02 Z©, the 
closest available to the mean metallicity in our simulations) as our 
default model (also see Figure 12). Then we run the MCRT code 
to compute the rate of ionizing photons that can escape the virial 
radius gesc- We define the escape fraction as /esc = 2esc/2int. 

4.1 Radiative Transfer Calculation 

We perform the MCRT using a Cartesian grid. We first convert 
each “well-resolved” galaxy identified in our simulations to a cu¬ 
bic Cartesian grid of side length L and with N cells along each 
dimension. We center the grid at the centre of the galactic halo 


and choose L equal to two virial radii. The size of a cell I — L/N 
must be appropriately chosen to ensure convergence. For each sim¬ 
ulation, we determine a minimum cell size /min and a maximum 
Amax and then take N — min{L//min, Anax}. We have run extensive 
tests to make sure the parameters /min and Anax for each simulation 
are carefully selected to ensure convergence for every snapshot and 
maintain reasonable computational expenses. We show examples 
of convergence tests in Appendix B^. These parameters are listed 
in Table 2. We then calculate the gas density, metallicity, and tem¬ 
perature, at each cell by distributing the mass, internal energy, and 
metals of every gas particle among a number of cells weighted by 
their SPH kernel. This conserves mass and energy of gas from the 
simulation to the grid. 

The MCRT method is similar to that described in Fumagalli et 
al. (2014). The radiation field is described by discrete Monte Carlo 
packets, each representing a large collection of photons of a given 
wavelength. We emit Astar packets isotropically from the location 
of the star particles, appropriately sampled by the star UV lumi¬ 
nosities. We also emit Auvb packets from the edge of the compu¬ 
tational domain in a manner that produces a uniform, isotropic UV 
background radiation field with intensity given by Faucher-Giguere 
et al. (2009). Every photon packet is propagated until it either es¬ 
capes the grid, or is absorbed somewhere in the grid. Scattering 
is included in the transport - i.e., we do not make the on the spot 
approximation. 

The photon packets are used to construct estimators of the hy¬ 
drogen photoionization rates in all cells. The photoionization cross- 
sections were taken from Verner et al. (1996), the collisional ion¬ 
ization rates from Jefferies (1968), and the radiative recombination 
rates from Verner & Eerland (1996). When calculating the rates, we 
use the gas temperature from the simulations instead of computing 
it self-consistently through the radiative transfer^. We use the case 
A recombination rates as the transport explicitly treats photon scat¬ 
tering. We assume that 40% of the metals are in dust phase and 
adopt a dust opacity of 10"^ cm^ g“^ (Dwek 1998; Eumagalli et al. 
2011). Since the high-redshift galaxies in our simulations tend to 
be extremely metal-poor, our results do not depend much on dust 
absorption. 

We assume that the gas is in ionization equilibrium, which 
should be valid for all but the lowest density, highest temperature 
regions. Such very low density regions likely do not influence the 
escape fraction in any case. We use an iterative method to reach 
equilibrium, running the MCRT, updating the ionization state of 
each cell, and then repeating the transport until convergence in the 
ionization state and escape fraction is reached. We use up to 15 

^ The MCRT calculation converges at much poorer resolution than that 
we use for hydrodynamics. This is because most of the sources reside in 
the environment where the ionizing photon optical depth is either tuv ^ 
1 or Tuv 1- The MCRT calculation will converge as long as the grid 
captures which limit a star particle is in. However, we emphasize that the 
high resolution of hydrodynamics is necessary in order to capture the ISM 
structure in star-forming regions in the presence of stellar feedback. Low 
resolution simulations tend to over-predict escape fractions by an order of 
magnitude (see the discussion in Section 5.2). 

^ The simulations take into account many other heating sources (e.g. 
shocks) besides photoionization. As the radiative transfer code includes col¬ 
lisional ionizations, it is more realistic to take the gas temperature from 
the simulations than re-computing gas temperature from radiative transfer 
calculations (in the latter case photoionization would be the only heating 
source). In regions dominated by photoionization, the uncertainty due to 
gas temperature is very small, since the recombination rate depends only 
weakly on temperature. 
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Figure 6. Stellar mass (top panels), star formation rate (second panels), escape fraction (third panels), and intrinsic and escaped ionizing photon budget 
(bottom panels, black and cyan lines, respectively) as a function of cosmic time for the most massive galaxy in each run. Open symbols show the time- 
averaged quantities over 100 Myr. The red dotted line in zSmlOmr shows the escape fraction calculated with the UV background turned off (Section 5.1). 
Instantaneous escape fractions are highly time variable, while the time-average escape fractions (over time-scales 100-1000 Myr) are modest (~ 5%). The 
intrinsic ionizing photon budgets are dominated by stellar population younger than 3 Myr, which tend to be embedded in dense birth clouds. Most of the 
escaping ionizing photons come from stellar populations aged 3-10 Myr, where a large fraction of sightlines have been cleared by stellar feedback. Note that 
the run using a common “sub-grid” star formation model (zSmlOe), which allows stars to form in diffuse gas instead of in dense clouds, severely over-estimates 

/esc- 
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log /esc 

Figure 7. Angular distribution of escape fraction for two typical snapshots, 
with spatially averaged escape fraction 0.005 (blue) and 0.2 (green), respec¬ 
tively. Statistics are obtained from N = 400 uniformly sampled directions. 
The broad distribution implies that the ionizing photons that escape to the 
IGM are highly anisotropic, and that the measured escape fraction from 
individual galaxies can vary by more than 2 dex depending on the sightline. 


iterations to reach convergence, with typical particle counts per it¬ 
eration of 3 X 10^ for A4tar and A^vb- We ran tests increasing the 
particle counts by an order of magnitude to check that the final es¬ 
cape fraction did not change. 
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4.2 Instantaneous and Time-averaged Escape Fraction 

In Figure 6, we present the instantaneous escape fractions (/esc), in¬ 
trinsic ionizing photon budgets (Qint), and escaped photon budgets 
(2esc) as a function of cosmic time for the most massive galaxies in 
each simulation. We also average Qm and Qesc over 100 Myr to ob¬ 
tain the time-averaged escape fractions (the open symbols in Figure 
6). The instantaneous escape fractions show significant time fiuctu- 
ations, varying between < 0.01% and > 20% from time to time. In 
our standard runs with default star formation prescriptions, galax¬ 
ies can reach high escape fractions (10-20%) only during small 
amounts of time. For most of the time, the time-averaged escape 
fractions remain below 5%. We also calculate the average escape 
fraction over their entire star formation history (i.e. z = 6-12). All 
our standard runs show values between 3-7%, which confirm low 
escape fractions on even longer time-scales. The variation in escape 
fractions on short time-scales is a consequence of feedback and the 
stochastic formation and disruption of individual star-forming re¬ 
gions, while long time-scale fluctuations are associated with galaxy 
mergers and intensive starbursts. Note that a high instantaneous es¬ 
cape fraction does not necessarily indicate a high contribution of 
ionizing photons. For example, although the main galaxy of z5mll 
had an escape fraction around 20% at z ~ 6.8, its intrinsic ionizing 
photon budget Qm was relatively low at that instant and the time- 
averaged escape fraction is only ~ 3%. Recalling that many models 
of reionization usually require /esc ~ 20% if the universe was reion¬ 
ized by galaxies brighter than Muv > —13 only (e.g. Finkelstein et 
al. 2012; Kuhlen & Faucher-Giguere 2012; Robertson et al. 2013), 
the escape fractions we find from our simulations are considerably 
lower than what those models require. 

As long as we properly resolve star formation in dense birth 
clouds, our results are not sensitive to the details of our star for¬ 


Figure 8. Time-averaged escape fraction (top panel) and escaped ionizing 
photon budget (bottom panel) as a function of cosmic time, color-coded by 
stellar mass. Different symbols represent the galaxies from different simu¬ 
lations. Points are the escape fraction averaged over 100 Myr (2esc/2int)- 
We see no significant dependence of /esc on redshift. 


mation prescription^. For example, in our zSmlOh run where we 
apply rtth = 1000 cm“^, the escape fraction is very similar to the 
standard runs. Also, the similarity between the HR z5ml0 run and 
the MR zSmlOmr run shows that our results converge with respect 
to resolution^. 

However, in our zSmlOe run where we allow stars form in dif¬ 
fuse gas, the time-averaged escape fraction exceeds 20% for most 
of the time. While this toy model results in higher escape fractions, 
we stress that such a star formation prescription is not consistent 
with our current understanding of star formation. As such, these 
predictions are likely not realistic but we include them to illus¬ 
trate how escape fraction predictions depend sensitively on the ISM 
model, with our zSmlOe run being representative of many simula- 


^ Previous studies also showed that GMC lifetimes and integrated star for¬ 
mation efficiencies were nearly independent of the instantaneous density 
threshold and star formation efficiency in dense gas, as long as the clouds 
were resolved (Hopkins et al. 2011, 2012). 

^ For HR runs where the mass of a star particle is only 10-100 Mq , we also 
test the effects of the IMF-average approximation. We randomly resample 
the ionizing photon budgets among individual star particles at a 1:20 ratio 
according to their ages and repeat the radiative transfer calculations. We 
find that the escape fractions are very similar. This confirms that the IMF- 
averaged approximation in HR runs does not affect our results. 
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Figure 9. Time-averaged escape fraction (top panel) and escaped ionizing 
photon budget (bottom panel) as a function of stellar mass, color-coded 
by cosmic time. Different symbols represent the galaxies from different 
simulations. Points are the escape fraction averaged over 100 Myr. The 
cyan dotted line in the bottom panel shows the best linear fit of log gesc = 
log(M*/M 0 ) +43.53. We see no strong dependence of /esc on M*. The 
dependence of 0esc on broadly follows the SFR-M* relation. 


tivities are intensive and can last for considerable amount of time, 
will the ionized regions expand and overlap and thus allow a large 
fraction of ionizing photons from the youngest stars to escape. For 
example, the high escape fractions in zSmlOmr at cosmic time > 1 
Gyr (z < 6 ) are due to the strong and lasting star formation during 
the past 100 Myr. However, such events are not common in our sim¬ 
ulations, since further star formation activity is usually suppressed 
effectively by stellar feedback. 

In Figure 7, we show the angular distribution of escape frac¬ 
tion as measured from N — 400 uniformly sampled directions. We 
repeat the radiative transfer calculation with ten times more photon 
packages than listed in Table 2 for two snapshots which have spa¬ 
tially averaged escape fraction /esc = 0.005 and 0.2, respectively. 
The broad distribution of escape factions implies that the ionizing 
photons escaping to the IGM from galaxies are highly anisotropic. 
It also indicates that the observationally measured escape fraction 
from individual galaxies does not necessarily reflect the angle aver¬ 
aged escape fraction from the same object, as it can vary by roughly 
2 dex depending on the sightline. 

In Figures 8 and 9, we compile the time-averaged escape 
fraction and escaped ionizing photon budget as a function of cos¬ 
mic time and stellar mass, respectively, for all the “well-resolved” 
galaxies in our standard runs. The symbols are color-coded by stel¬ 
lar mass and cosmic time in Figure 8 and 9, respectively. Most of 
the points lie below /esc <5%. The escape fraction has a large scat¬ 
ter at flxed cosmic time or stellar mass. We And that there is no 
signiflcant dependence of escape fraction on cosmic time or stellar 
mass. This is consistent with the argument that the escape of ion¬ 
izing photons is restricted by small-scale ISM structures surround¬ 
ing the young stellar populations. More simulations are required to 
study possible redshift and galaxy mass evolution to lower redshifts 
and over a wider mass interval than sampled by the simulations an¬ 
alyzed in this paper. We also caution that weak trends would be dif- 
flcult to discern given the time variability found in our simulations 
and the small size of our simulation sample. The escaped ionizing 
photon budget depends linearly on stellar mass, with the best fit 
logfiesc = log(^*/M 0 ) +43.53. This is primarily a consequence 
of the roughly linear dependence of SFR on stellar mass. 


tions that do not have sufficient resolution to capture dense ISM 
structures. 

For the flducial stellar population model we adopt, the major¬ 
ity of the intrinsic ionizing photons are produced by the youngest 
star particles with age < 3 Myr (see also Figure 12). These stars are 
formed in dense, self-gravitating, molecular regions. Most of their 
ionizing photons are immediately absorbed by their “birth clouds” 
and thus cannot escape the star-forming regions (see also, e.g. Kim 
et al. 2013). When a star particle is older than 3 Myr, a large cov¬ 
ering fraction of its birth cloud has been cleared by feedback and 
thus a signiflcant fraction (order unity) of its ionizing photons are 
able to propagate to large distances (see e.g. the middle panels of 
Figure 2). Indeed, the ionizing photons that escaped in our simu¬ 
lations mostly come from the star particles of age between 3-10 
Myr (also see Section 5.4). However, the intrinsic ionizing photon 
budget of a star particle decreases rapidly with age above 3 Myr ac¬ 
cording to many standard stellar population models. In other words, 
the escape fractions are primarily determined by small-scale ISM 
structures surrounding young and intermediate-age star particles. 
The low escape fractions we And in our simulations are the conse¬ 
quence of the fact that the time-scale for a star particle to clear its 
birth cloud is comparable to the time-scale for it to exhaust a large 
amount of its ionizing photon budget. Only when star formation ac- 


5 DISCUSSION 

We And that instantaneous escape fractions of hydrogen ionizing 
photons from our simulated galaxies vary between 0 . 01 %- 20 % 
from time to time, while time-averaged escape fractions generally 
remain below 5%. These numbers are broadly consistent with the 
wide range of observationally constrained escape fractions mea¬ 
sured from variant galaxy samples at z = 0-3 (e.g. Leitet et al. 2011, 
2013; Cowie et al. 2009; Siana et al. 2010; Bridge et al. 2010; Iwata 
et al. 2009; Boutsia et al. 2011; Vanzella et al. 2012; Nestor et al. 
2013). 

We obtain much lower escape fractions than previous simu¬ 
lations with “sub-grid” ISM, star formation, and feedback models 
(e.g. Razoumov & Sommer-Larsen 2010; Yajima et al. 2011), but 
our results are more consistent with many recent simulations with 
state-of-art ISM and feedback models (e.g. Kim et al. 2013; Wise 
et al. 2014; Kimm & Cen 2014; Paardekooper et al. 2011, 2015). 
Below, we will show that this owes to the failure of the “sub-grid” 
models in resolving stellar birth clouds. 

Nevertheless, the escape fractions from our simulated galax¬ 
ies are still considerably lower than what requires for these galaxies 
to reionize the universe in many popular models. The tension can 


© 0000 RAS, MNRAS 000, 1-17 








12 Ma et al. 


be at least partly resolved by invoking galaxies much fainter than 
what we study in this work, since smaller galaxies have dramati¬ 
cally increasing number densities and possibly much higher escape 
fractions (e.g. Alvarez et al. 2012; Paardekooper et al. 2013; Wise 
et al. 2014). Alternatively, in the rest of this section, we will discuss 
some physical parameters that might boost the escape fractions in 
our simulated galaxies. Most of the experiments presented here are 
for illustrative purposes, but they are worth further exploration in 
future work in a more systematic and self-consistent way. 

5.1 UV Background 

We repeat the radiative transfer calculation for all the snapshots af¬ 
ter cosmic time 0.9 Gyr (z ~ 6) of our zSmlOmr galaxy with the 
UV background switched off (the red dotted line in the upper right 
panel of Figure 6). The predicted escape fractions does not differ 
from the previous calculation with the UV background at 0.01% 
level, consistent with the results in Yajima et al. (2011). This con¬ 
firms that the low-density, diffused gas in the galactic halo (which 
is affected by the UV background) does not affect much the escape 
of ionizing photons 

5.2 Star Formation Criteria 

In our standard simulations, we allow star formation occurs only in 
molecular, self-gravitating gas with density above a threshold Uth = 
100 cm“^. We run z5ml0h where we adopt = 1000 cm“^ for 
a convergence study. For contrast, we intentionally design z5ml0e 
to mimic “sub-grid” SF models, where we lower nth to 1 cm“^ and 
allow extra SF at 2% efficiency per free-fall time in gas above the 
threshold but not self-gravitating. In Section 3, we have confirmed 
that the global galaxy properties (e.g. star formation rates, stellar 
masses, UV magnitudes, etc.) are very similar between these runs. 
However, as it shown in Figure 6, the escape fraction from z5ml0e 
is significantly higher than in other simulations. 

To illustrate this more clearly, we compare the time-averaged 
(over 100 Myr time-scale) escape fraction of z5ml0, z5ml0e, and 
z5ml0h in Figure 10. The qualitative behaviors of escape fraction 
are very similar between z5ml0, z5ml0mr, and z5ml0h, which 
further confirm that our results are converged with respect to res¬ 
olution and SF density threshold (as long as it is much larger than 
the mean density of the ISM). 

However, in z5ml0e, the escape fraction is dramatically 
higher, since many young stars form in the diffuse ISM. Their ion¬ 
izing photons can then immediately escape the galaxy. We empha¬ 
size that the z5ml0e run is not realistic but mimics “sub-grid” SF 
models as adopted in low-resolution simulations where star forma¬ 
tion in dense gas clouds cannot be resolved. This suggests a caution 
that simulations with “sub-grid” SF models can over-predict the es¬ 
cape fraction by an order of magnitude. 

5.3 Runaway Stars 

There is plenty of evidence that a considerable fraction of O and B 
stars have high velocities and can travel far from their birth clouds 
during their lifetime (the “runaway” stars, e.g. Blaauw 1961; Stone 
1991; Hoogerwerf et al. 2001; Tetzlaff et al. 2011). To qualitatively 
illustrate the effect of these runaway stars on the escape fraction 

However, if the simulations are run without a UV background, gas ac¬ 
cretion onto the halo itself can be modified. 
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Figure 10. Escape fraction with different star formation density prescrip¬ 
tions. The escape fractions averaged over 100 Myr are shown for z5ml0 
(72th = 100 cm“^, blue solid), zSmlOmr ( 22 th = 100 cm“^, green dotted), 
zSmlOh (72th = 1000 cm“^, cyan dashed), and zSmlOe ( 72 th = 1 cm“^, 
without SF self-gravity criteria, red dotted). For 72th ^100 cm“^, our re¬ 
sults are well-converged with respect to SF density threshold and resolution. 
However, if SF is allowed in diffuse gas, /esc can be severely over-estimated. 

(e.g. Conroy & Kratter 2012), we move every star particle younger 
than 3 Myr by a distance Vini4ge along a random direction in the 
snapshots, and repeat the radiative transfer calculation to evaluate 
the escape fraction as the stars are at their new positions. Here Vim is 
some initial kicking velocity and tage is the age of the star particle. 
In principle, it would be more self-consistent if we re-run the whole 
simulation with runaway star prescription (e.g. Kimm & Cen 2014). 
Nonetheless, our simple experiment provides a first estimate of the 
effects of runaway stars on the escape fraction. 

We repeat this experiment for our z5ml0mr run during cos¬ 
mic time between 0.8-1.0 Gyr (z = 6-7) with vm = 10, 50, and 
100 km s“\ corresponding to a displacement of 30, 150, and 300 
pc for a star particle of age 3 Myr. We show the results in Figure 
11. A small initial velocity of Vmi = 10 km s“^ barely affects the es¬ 
cape fraction since the displacement of a newly formed star particle 
is < 30 pc, which is much less than the typical size of their birth 
clouds (see Figure 2 for an illustration of the ISM structure around 
young star particles). For vm = 50 km s“\ the escape fractions can 
be boosted by at most 1-2% (in absolute units, or 20-30% fraction¬ 
ally). Only for extremely high initial velocity (~ 100 km s“^), the 
escape fractions are enhanced by a few precent, since some young 
star particles are kicked out of their birth clouds. But these num¬ 
bers are still somewhat lower than what many reionization models 
require. Note that observations suggest that only ~ 30% of the OB 
stars are runaway stars and that the typical velocity of runaway stars 
is around 30 km s~^ (e.g. Tetzlaff et al. 2011). Therefore, our ex¬ 
periments suggest that runaway stars will boost the escape fractions 
by no more than 1% (in absolute units, or 20% fractionally) 

One effect not captured in our post-processing experiment which could 
potentially boost the escape fraction is how feedback from runaway stars 
would affect the structure of the ISM as they move away from their birth 
locations. A self-consistent modeling of runaway stars is presented in Kimm 
& Cen (2014). Our result is broadly consistent with theirs for halo masses 
- 10^0 M©. 
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Figure 11. Escape fractions in presence of runaway stars. We only show 
zSmlOmr during the cosmic time 0.8-1.0 Gyr (z = 6-7, but the effect 
in other runs would is similar). Each star particle younger than 10 Myr 
is kicked from its original position along a random direction with an ini¬ 
tial velocity Vini. Blue dotted, dashed, and solid lines show the results for 
Vjni = 10kms“^, 50kms“^, and 100kms“^, respectively. The black 
solid line shows the escape fraction when vini = 0 (the same as in Fig¬ 
ure 6). Typical kick velocities suggested by observations (~ 30 km s“^) 
have only small effects on /esc- Only if the velocities are very large (e.g. 
>100 km s“^), and an order-unity fraction of stars have been kicked, will 
this be significant. 


5.4 Stellar Population Models 

So far, we have adopted the Padova+AGB stars track of metallic- 
ity Z = 0.0004 (0.02 Z©) from STARBURST99 model assuming a 
Kroupa (2002) IMF from 0.1-100 M© to evaluate the intrinsic ion¬ 
izing photon budget for each star particle. In this model, the ion¬ 
izing photon production rate decreases rapidly when the age of a 
stellar population exceeds 3 Myr. However, there are good reasons 
to believe that these models suffer from great uncertainties. For ex¬ 
ample, Steidel et al. (2014) emphasized the importance of binary 
and rotating stars since these stars have high effective temperatures 
that are required to explain the ionization states of z = 2-3 star 
forming galaxies. Moreover, recent theoretical studies suggest that 
binary star interactions can produce ionizing photons in a stellar 
population older than 3 Myr; such events are not uncommon (e.g. 
de Mink et al. 2014). While these models are very uncertain and 
still poorly understood, it is not unphysical to invoke more ionizing 
photons from these populations. To explore the effects of different 
stellar population models on the escape fractions, we construct a 
toy model which we refer to the “constant gint model” to explore 
its effect on the escape of ionizing photons. In this model, the ion¬ 
izing photon budget of a stellar population is 5.6x Mo“‘ 

when the population is younger than 10 Myr and suddenly drops 
to 5.6 X 10^^ s~^ M©^^ when the population is older than 10 Myr. 
For comparison, we also tabulate the ionizing photon budget using 
the Padova-fAGB stars track of solar metallicity (Z = 0.02) from 
STARBURST99 model. We illustrate in Figure 12 the intrinsic ion¬ 
izing photon budget as a function of stellar age for the three models 
we discuss. Their behaviors are very similar for stellar age < 3 Myr, 
after which they start to deviate heavily from each other. The solar- 
metallicity model has the lowest ionizing photon budget between 
3-10 Myr while the constant gmt model has the highest. 



Age (Myr) 


Figure 12. Ionizing photon budget per unit mass for a stellar population 
as a function of its age. The black line shows the STARBURST99 low- 
metallicity model (our default model). The green line shows the STAR- 
BURST99 solar-metallicity model. The red line shows a simple “constant 
0int model” we consider in Figure 13. This produces a similar number of 
ionizing photons in the first 3 Myr, but retains the same photon production 
rate to 10 Myr. 



Cosmic time (Gyr) 

Figure 13. Escape fractions calculated by invoking the three different stellar 
population models shown in Figure 12. We show the results for zSmlOmr 
during cosmic time 0.8-1.0 Gyr (z = 6-7, but the effect in all other runs is 
similar). The black line shows the results using the STARBURST99 low- 
metallicity model (Z = 0.0004, our fiducial model, the same as in Figure 6). 
The green line shows the results using STARBURST99 solar-metallicity 
model (Z = 0.02). The red line shows the results when using the “constant 
0int model”. By extending the lifetime of photon production to 3-10 Myr, 
when the birth clouds have been largely cleared, large /esc (10-20%) can be 
obtained. 

We repeat the radiative transfer calculation to calculate the es¬ 
cape fraction assuming intrinsic ionizing photon production rate 
evaluated from these models. In Figure 13, we show the results for 
our zSmlOmr run during cosmic time 0.8-1.0 Gyr (z = 6-7). The 
escape fractions are very sensitive to the stellar models we use. For 
the solar-metallicity track model, we get the lowest escape frac- 
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tions. On the other hand, if we adopt the constant Qmi model, we 
find the escape fractions are enhanced by almost an order of mag¬ 
nitude. These results further illustrate the picture that the escaped 
photons come from star particles of age 3-10 Myr, where their 
birth clouds have been cleared by feedback. Our findings suggest 
that relatively older stellar populations could contribute a consider¬ 
able fraction of ionizing photons during reionization, if these pop¬ 
ulations produce more ionizing photons than what standard stellar 
population models predict, as motivated by models that include ro¬ 
tation, binaries, and mergers. 


6 CONCLUSIONS 

In this work, we present a series of extremely high-resolution (par¬ 
ticle mass 2O-2OOOM0, smoothing length 0.1-4 pc) cosmological 
zoom-in simulations of galaxy formation down to z ~ 6, covering 
galaxy halo masses in 10^-10^^ M©, stellar masses in 2 x 10^- 
2 X 10^ M©, and rest-frame ultraviolet magnitudes Muv = — 9 to 
— 19 at that time. This set of simulations include realistic models 
of the multi-phase ISM, star formation, and stellar feedback (with 
no tuned parameters), which allow us to explicitly resolve small- 
scale ISM structures. Cosmological simulations with these feed¬ 
back models have been shown to produce reasonable star forma¬ 
tion histories, the stellar-mass halo mass relation, the Kennicutt- 
Schmidt law, the star-forming main sequence, etc., at z = 0-6 (Hop¬ 
kins et al. 2014). We post-process our simulations with a Monte 
Carlo radiative transfer code to evaluate the escape fraction of hy¬ 
drogen ionizing photons from these galaxies. Our main conclusions 
include the following. 

(i) Instantaneous escape fractions have large time variabilities, 
fiuctuating from < 0.01% to > 20% from time to time. In our stan¬ 
dard runs, the escape fractions can reach 10-20% only for a small 
amount of time. The time-averaged escape fractions (over time- 
scales 100-1000 Myr) generally remain below 5%, considerably 
lower than many recent models of reionization require. 

(ii) As long as star formation is regulated effectively via feed¬ 
back, the escape fractions are mainly determined by small-scale 
ISM structures around young and intermediate-age stellar popula¬ 
tions. According to standard stellar population models, most of the 
intrinsic ionizing photons are produced by newly formed star parti¬ 
cles younger than 3 Myr. They tend to be embedded in their dense 
birth clouds, which prevent nearly all of their ionizing photons 
from escaping. The escaping ionizing photons primarily come from 
intermediate-age stellar populations between 3-10 Myr, where the 
dense birth clouds have been largely destroyed by feedback. Ac¬ 
cording to “standard” stellar population models, the ionizing pho¬ 
ton production rates decline heavily with time at these ages. This 
leads to the difficulty of getting high escape fractions. 

(iii) The escape fractions do not change if the star formation 
density threshold increases from 100 to 1000 cm“^, as long as stars 
form in resolved, self-gravitating, dense clouds. On the other hand, 
if we allow star formation in the diffuse ISM (with some ad hoc low 
star formation efficiency), as is adopted in most low-resolution cos¬ 
mological simulations, the escape fractions can be over-predicted 
by an order of magnitude. We emphasize that realistic, resolved 
phase structure of the ISM is critical for converged predictions of 
escape fractions. 

(iv) Applying a fraction of ~ 30% runaway OB stars to our 
simulations with typical velocity ~ 30 kms“^ as motivated by 
many observations can only enhance the escape fraction by at most 
1% (in absolute values, or 20% fractionally). The effect of run¬ 


away stars would not be significant unless a large fraction of the 
most young stars can obtain dramatically high initial velocity in 
high-redshift galaxies. 

(v) Stellar populations older than 3 Myr may play an impor¬ 
tant role in reionizing the universe. The escape fractions can be 
boosted significantly if stellar populations of intermediate ages pro¬ 
duce more ionizing photons than what standard stellar population 
models predict, as suggested by many new stellar population mod¬ 
els (e.g. models including rotations, binary interactions, and merg¬ 
ers). 

Our simulations are limited in sample size. Also, the sim¬ 
ple experiments we present in Section 5 do not treat stellar feed¬ 
back consistently with varying stellar models. Our results motivate 
further work exploring the effects of IMF variations, stellar evo¬ 
lution models, runaway stars, etc., in a more systematic and self- 
consistent way. 
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APPENDIX A: COMPARISON OF MCRT 

POST PROCESSING TO ON-THE-FLY IONIZATION 

CALCULATIONS 

As described in Section 2, in our simulations, the ionization state of 
each gas particle at every timestep is determined from the photoion¬ 
ization equilibrium equations described in Katz et al. (1996), given 
a uniform and redshift-dependent UV baekground from Eaucher- 
Giguere et al. (2009) and photon-ionization and photon-eleetric 
heating rate from local sources, assuming a loeal Jeans-length ap¬ 
proximation of self-shielding. In the simulations (“on-the-fly”), we 
model photoionization feedback from star particles in an approx¬ 
imate way - we move outward from the star particle and ionize 
each nearest neutral gas particle until the photon budget is com- 
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Figure 1. Gas neutral fraction on a slice across the galactic centre. We show the snapshot at z = 6 from zSmlOmr. Left: the neutral fraction directly extracted 
from the simulation. Right: the neutral fraction recomputed by radiative transfer code. The size of the box equals to two times of the virial radius and there are 
250 pixels along each direction. Only regions within the virial radius are shown. 
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Figure 2. Neutral hydrogen column density as viewed from the galactic centre. We show the snapshot at z = 6 from zSmlOmr. Left: the neutral column density 
directly extracted from the simulation. Right: the neutral column density recomputed by radiative transfer code. 



Cell size (pc) 

Figure 3. Resolution convergence of MCRT calculation. We show two ex¬ 
amples, one from z5m09 (HR) and the other from zSmlOmr (MR), where 
we repeat the MCRT calculation for the same galaxy but with different grid 
cell sizes for the RT calculation. The results are relatively insensitive to the 
cell size, for cell sizes varying from I = 20-100 pc. 


pletely consumed. In intense star-forming regions, this allows H II 
regions to expand and overlap and thus approximately captures rea¬ 
sonable ionization states in these regions. However, if the gas dis¬ 
tribution is highly asymmetric around an isolated star particle (see, 
e.g. the middle column in Figure 2), the gas ionization states will 
not be accurately captured. In this work, we follow the propagation 
of ionizing photons and re-compute the gas ionization state with a 
Monte Carlo radiative transfer code in post-processing, which will 
be more accurate in photoionization regions. 

Figure 1 shows the gas neutral fractions on a slice crossing 
the galactic centre, of a snapshot at z = 6 from z5ml0. Figure 2 
shows the neutral hydrogen column density map as viewed from 
the centre of the galaxy for the same snapshot. In both figures, the 
left panels are the results before post-processing and the right pan¬ 
els show the results from radiative transfer calculations (using ten 
times the standard number of photon packages listed in Table 2). 
In general, both results agree quite well on large-scale pattern of 
the neutral gas distribution, although radiative transfer calculations 
reveal more small structures in star-forming regions. None of our 
conclusions in this paper are changed qualitatively if we compute 
the escape fractions using on-the-fiy ionization states in the simu¬ 
lations. It is reassuring, both for the present work and for previous 
studies that used the same approximations, that the approximations 


© 0000 RAS, MNRAS 000, 1-17 















Escape Fractions from FIRE Galaxies 


17 


used in the simulation code predict ionization structures that are 
broadly consistent with post-processing radiative transfer calcula¬ 
tions. 


APPENDIX B: RESOLUTION CONVERGENCE FOR 
MCRT CALCULATION 

In Section 4.1, we describe that the MCRT calculation is performed 
on a cubic Cartesian grid of side length L and with N cells along 
each dimension. In principle, the resolution I — L/N should be 
small enough to capture the ISM structure, but the number of cells 
cannot be so big that the calculation is too computationally ex¬ 
pensive. After performing extensive convergence tests, we choose 
/ = 25-100 pc depending on the size of the galactic halo. Here we 
show two typical examples, one from z5m09 (HR) and the other 
from zSmlOmr (MR), to illustrate the convergence of our MCRT 
calculation with respect to cell size /. As shown in Figure 3, we 
repeat the MCRT calculation for the same galaxy with resolution 
varying from / = 20-100 pc and find that the escape fractions do 
not change appreciably. 

The MCRT calculation converges at much poorer resolution 
than that we adopt for hydrodynamics. This is because most of the 
sources reside in an environment where the ionizing photon opti¬ 
cal depth is either tuv ^ 1 or tuv <C 1. In both limits, the MCRT 
calculation converges even if the exact column density is not cap¬ 
tured with great accuracy (e.g. tuv = 100 and tuv = 10 make lit¬ 
tle difference). However, we emphasize that the high resolution of 
hydrodynamics is necessary in order to capture the ISM structure 
in the presence of star formation and stellar feedback. Low resolu¬ 
tion simulations with “sub-grid” models tend to over-predict escape 
fraction (see the discussion in Section 5.2). 
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